#pragma once

#include "matrix.h"




template <int TEMPL_DEGREE, class ARG_TYPE>
class polynomial
{
	ARG_TYPE m_fScale;
	ARG_TYPE m_fScaleInv;
	mdsVector<TEMPL_DEGREE, ARG_TYPE> m_vecCoeffs;

public:

	polynomial() { m_fScaleInv = m_fScale = 1; }

	polynomial(const polynomial<TEMPL_DEGREE, ARG_TYPE> & polySrc )
	{
		m_vecCoeffs = polySrc.m_vecCoeffs;
		m_fScaleInv = polySrc.m_fScaleInv; m_fScale = polySrc.m_fScale;
	}

	operator=(const polynomial<TEMPL_DEGREE, ARG_TYPE> & polySrc )
	{
		m_vecCoeffs = polySrc.m_vecCoeffs;
		m_fScaleInv = polySrc.m_fScaleInv; m_fScale = polySrc.m_fScale;
	}

	const mdsVector<TEMPL_DEGREE, ARG_TYPE> & GetCoeffVec() const { return m_vecCoeffs;}

	void SetCoeffVec(const mdsVector<TEMPL_DEGREE, ARG_TYPE> & vecSetTo)
	{ 
		m_vecCoeffs = vecSetTo;
	}

	ARG_TYPE & Coeff(int iWhichCoeff)
	{
		return m_vecCoeffs[iWhichCoeff];
	}

	const ARG_TYPE & Coeff(int iWhichCoeff) const
	{
		return m_vecCoeffs[iWhichCoeff];
	}

	void Differentiate() // works in-situ, clobbers old values
	{
		for(int i = 1; i < TEMPL_DEGREE; ++i)
		{
			m_vecCoeffs[i-1] = i*m_vecCoeffs[i];
		}
		m_vecCoeffs[TEMPL_DEGREE] = 0;
	}

	void ScaleWidthBy(const ARG_TYPE & fScaleBy)
	{
		ARG_TYPE fScaleInv = 1/fScaleBy;
		m_fScale *= fScaleBy; m_fScaleInv *= fScaleInv;
		ARG_TYPE fCoeffScale = 1;
		for(int i = 0; i < TEMPL_DEGREE; ++i)
		{
			m_vecCoeffs[i] *= fCoeffScale;
			fCoeffScale *= fScaleInv;
		}
	}

	void ShiftBy(const ARG_TYPE & fShiftBy)
	{
		ARG_TYPE arrFPascalCoeffs[TEMPL_DEGREE];
		for(int i = 0; i < TEMPL_DEGREE; ++i)
		{
			arrFPascalCoeffs[i] = 0;
		}
		int nIdxLevel;
		for(nIdxLevel = 0; nIdxLevel < TEMPL_DEGREE; ++nIdxLevel)
		{

			// add shift to smaller-degree coeffs
			ARG_TYPE fShiftPower = -fShiftBy;
			for(int iPower = nIdxLevel-1; iPower >= 0; --iPower)
			{
				m_vecCoeffs[iPower] =
					m_vecCoeffs[iPower]+
					fShiftPower*arrFPascalCoeffs[nIdxLevel-iPower-1]*m_vecCoeffs[nIdxLevel];
				fShiftPower *= -fShiftBy;
			}
			
			// update pascal's triangle coefficients
			for(int nIdxCoeff = nIdxLevel; nIdxCoeff >0 ; --nIdxCoeff)
			{
				arrFPascalCoeffs[nIdxCoeff] = arrFPascalCoeffs[nIdxCoeff]+
					arrFPascalCoeffs[nIdxCoeff-1];
			}
			arrFPascalCoeffs[0] += 1;
		}
	}

	ARG_TYPE GetScale() const { return m_fScale; }
	ARG_TYPE GetScaleInv() const { return m_fScaleInv; }

	void ScaleHeightBy(const ARG_TYPE & fScaleHeight)
	{
		for(int i = 0; i < TEMPL_DEGREE; ++i)
		{
			m_vecCoeffs[i] *= fScaleHeight;
		}
	}

	ARG_TYPE Integrate(const ARG_TYPE & fConstant) // WARNING: read other comments thoroughly
	{
		// works in-situ, clobbers old values, return value is highest-degree term
		// note that integration adds a degree to a polynomial.  This run-over
		// is stored in the return value, but not internally in the polynomial.
		// If highest-degree coefficient is zero, integration is safe.
		ARG_TYPE fResult = m_vecCoeffs[TEMPL_DEGREE-1]/TEMPL_DEGREE;
		for(int i = TEMPL_DEGREE-1; i > 0; --i)
		{
			m_vecCoeffs[i] = m_vecCoeffs[i-1]/i;
		}
		
		m_vecCoeffs[0] = fConstant;
		return fResult;
	}

	ARG_TYPE IntegrateIval(const ARG_TYPE & fMin, const ARG_TYPE & fMax) const
	{
		ARG_TYPE fResultMin = 0.0;
		
		for(int i = TEMPL_DEGREE-1; i >= 0; --i)
		{
			fResultMin = fMin*(fResultMin+m_vecCoeffs[i]/(i+1.0));
		}

		ARG_TYPE fResultMax = 0.0;
		for(int i = TEMPL_DEGREE-1; i >= 0; --i)
		{
			fResultMax = fMax*(fResultMax+m_vecCoeffs[i]/(i+1.0));
		}
		return fResultMax-fResultMin;
	}

	ARG_TYPE EvalUnscaled(const ARG_TYPE & fTin) const
	{
		ARG_TYPE fT = fTin * m_fScale;
		ARG_TYPE fTPow = 1;
		ARG_TYPE fAccum = 0;
		for(int i = 0; i < TEMPL_DEGREE; ++i)
		{
			fAccum += fTPow*m_vecCoeffs[i];
			fTPow *= fT;
		}
		return fAccum;
	}

	ARG_TYPE EvalScaled(const ARG_TYPE & fT) const
	{
		ARG_TYPE fTPow = 1;
		ARG_TYPE fAccum = 0;
		for(int i = 0; i < TEMPL_DEGREE; ++i)
		{
			fAccum += fTPow*m_vecCoeffs[i];
			fTPow *= fT;
		}
		return fAccum;
	}
		
};

namespace polynomialHelpers
{
	template <int TEMPL_DEGREE_OUT, int TEMPL_DEGREE_IN, class ARG_TYPE>
	polynomial<TEMPL_DEGREE_OUT, ARG_TYPE> promoteUp(const polynomial<TEMPL_DEGREE_IN, ARG_TYPE> & src)
	{
		polynomial<TEMPL_DEGREE_OUT, ARG_TYPE> result;
		if(TEMPL_DEGREE > TEMPL_DEGREE_OUT)
		{
			ASSERT(0); // I wish I knew of a way to check this at compile time
		}
		int i = 0;
		for(; i < TEMPL_DEGREE_IN; ++i)
		{
			result.Coeff(i) = src.Coeff(i);
		}
		for(; i < TEMPL_DEGREE_OUT; ++i)
		{
			result.Coeff(i) = static_cast<ARG_TYPE>(0);
		}
		return result;
	}

	template<int TEMPL_DEG_IN1, int TEMPL_DEG_IN2, class ARG_TYPE>
	polynomial<TEMPL_DEG_IN1+TEMPL_DEG_IN2, ARG_TYPE> multiply(
		const polynomial<TEMPL_DEG_IN1, ARG_TYPE> & pIn1,
		const polynomial<TEMPL_DEG_IN2, ARG_TYPE> & pIn2)
	{
		polynomial<TEMPL_DEG_IN1+TEMPL_DEG_IN2, ARG_TYPE> result;
		for(int i = 0; i < TEMPL_DEG_IN1; ++i)
		{
			for(int j = 0; j < TEMPL_DEG_IN2; ++j)
			{
				result.Coeff(i+j) = result.Coeff(i+j)+
					p1.Coeff(i)*p2.Coeff(j);
			}
		}
		return result;
	}

};

#define MAX_SPLINE_POLY_DEGREE 4

typedef polynomial<MAX_SPLINE_POLY_DEGREE, double> SplineEditCurve;

